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This  memorandum  describes  progress  in  the  development  of  a set  of  computer 
programs  which  solve  the  equations  of  gas/particle  flow  through  an  axisymmetric 
nozzle.  Two  programs  have  been  written  to  solve  the  equations  in  the  transonic 
throat  region  of  the  nozzle.  The  first  treats  the  two-phase  fluid  as  a heavy 
perfect  gas  with  modified  isentropic  exponent  and  molecular  weight,  and  solves 
the  equations  of  isentropic  transonic  flow.  An  initial  flow  field  configuration 
is  obtained  for  input  to  the  second  program,  which  incorporates  non-equilibrium 
effects  of  the  gas/particle  mixture. 

The  programs  provide  accurate  initial  line  data  for  starting  a supersonic 
calculation,  and  will  be  used  as  input  to  a third  program,  soon  to  be  written, 
for  solving  the  two-phase  flow  equations  in  the  supersonic  region  of  the  nozzle. 
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1 


INTRODUCTION 


Accurate  prediction  of  the  properties  of  a rocket  exhaust  plume  necessi- 
tates detailed  analyses  of  the  nozzle  and  plume  flow  fields.  The  Rocket  Exhaust 
Plume  program  (REP)1  and  the  Base  Flow  program  (BAFL)2,  which  calculate  the 
structure  of  the  exterral  (plume)  flow  field,  incorporate  the  latest  develop- 
ments in  gas  dynamics  and  chemical  kinetics,  and  are  the  most  advanced  models 
available.  However,  the  accuracy  of  the  predictions  depends  upon  the  accuracy 
of  the  input  nozzle  exit  parameters,  and  the  programs  currently  used  at  the  RPE 
for  calculating  the  flow  from  the  combustion  chamber  to  the  nozzle  exit  plane3, 
while  able  to  take  account  of  non-equilibrium  chemical  reactions,  can  treat  only 
one-dimensional,  gas  phase  flow. 

Many  investigators  (e.g.  Ref.  5)  have  shown  that  the  flow  in  a typical 
rocket  nozzle  is  far  from  one-dimensional,  and  it  has  also  been  shown  (e.g. 

Ref.  6)  that  condensed  metal  oxides  in  the  rocket  exhaust  nozzle  can  signifi- 
cantly affect  the  performance  of  the  motor.  These  factors  must  be  considered 
when  the  nozzle  calculations  are  performed,  and  the  ultimate  aim  of  current  work 
is  to  develop  a set  of  computer  programs  to  solve  the  equations  governing  the 
flow  of  a gas/particle  mixture  through  an  axisymmetric  rocket  exhaust  nozzle. 

In  the  solving  of  the  equations  for  flow  through  a convergent/divergent 
nozzle,  many  mathematical  difficulties  are  encountered  in  the  region  of  transi- 
tion from  subsonic  to  supersonic  flow.  Although  the  transition  is  smooth  the 
properties  in  the  two  regions  are  quite  different.  In  the  subsonic  region  the 
governing  equations  are  elliptic,  whereas  in  the  supersonic  region  they  are 
hyperbolic.  The  supersonic  flow  field  may  be  solved  by  the  method  of  character- 
istics, but  accurate  initial  line  data  are  needed  to  start  the  calculation. 

Since  the  flow  near  the  nozzle  throat  is  far  from  one -dimensional , the  sonic 
velocity  at  the  wall  being  reached  well  upstream  of  the  throat,  while  that  at 
the  axis  is  reached  downstream  of  the  throat,  a method  of  solution  which 
describes  the  properties  of  both  subsonic  and  supersonic  flow  segments  must  be 
devised.  This  memorandum  describes  the  method  adopted  for  solving  the  flow  in 
this  transonic  throat  region. 

2 TRANSONIC  TWO-PHASE  FLOW  IN  AXISYMMETRIC  NOZZLES 

"To  solve  the  equations  governing  the  transonic  flow  of  a gas-particle 
mixture  is  an  extremely  formidable  task"  - Kliegel6. 

Of  the  many  approaches  which  have  been  made  to  the  solution  of  the  two- 
phase  flow  problem1*,  the  earliest  was  to  treat  the  nozzle  expansion  processes  as 

. - - • ■ • 1 — '■  1 ■ — i 

PFECEDIMJ  PAOE|)BLANK-NOT  FILMED 


5 


uncoupled  - that  is,  the  particle  velocity  and  thermal  lags  were  considered  to 
be  independent  of  each  other,  and  their  effects  on  the  gas  properties  were 
ignored.  More  recently,  Kliegel  and  Nickerson7  have  developed  the  method  of 
constant  fractional  lag,  in  which  the  ratio  of  the  gas  to  particle  velocities 
and  the  ratio  of  the  temperatures  are  represented  by  constants.  Their  program 
is  the  most  widely  used,  but  the  method  of  solution  is  restricted  to  nozzles  of 
a particular  shape,  and  is  inapplicable  to  nozzles  of  sharp  wall  curvature. 

It  was  decided  to  use  the  technique  employed  by  Regan,  Thompson  and 
Hoglund9,  as  this  leads  to  a fully  coupled  solution,  and  is  valid  for  a much 
wider  range  of  nozzle  shapes.  The  method  proceeds  in  two  stages,  and  a separate 
computer  program  has  been  written  for  each. 

In  the  first  stage,  the  gas/particle  mixture  is  assumed  to  be  in  equili- 
brium - that  is,  there  are  no  particle  heat  or  velocity  lags,  and  the  two-phase 
fluid  is  treated  as  a heavy  perfect  gas  with  modified  molecular  weight 

m * m (1  + <J>)  (2.1) 

based  on  the  particle-to-gas  mass  ratio  <)>  , and  effective  isentropic  exponent 

Y * Y (1  + ♦ )/( 1 + Y<J>  ) (2.2) 

c c 


where 


0>  cPp/cPg 


The  equations  for  isentropic  transonic  flow  are  solved,  and  initial  flow 
field  data  obtained  for  input  to  the  second  stage.  The  equilibrium  program  may 
be  used  separately  to  provide  starting  values  for  a supersonic  solution  in  the 
absence  of  particles,  or  when  the  non -equilibrium  particle  effects  can  be 
ignored. 

In  the  second  program,  the  non-equilibrium  effects  of  the  two-phase  fluid 
are  incorporated.  The  partial  differential  equations  governing  the  gas/particle 
flow  are  rewritten  as  algebraic  replacement  equations  and,  using  the  results 
from  the  initial  stage  as  a first  approximation,  are  solved  iteratively. 
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STAGE  1:  TRANSONIC  EQUILIBRIUM  PARTICLE/GAS  MIXTURE  PROGRAM  (TEPGAS) 


There  are  several  recognised  methods  for  solving  the  equations  of  transonic, 
axisymmetric  flow,  but  most  are  applicable  only  to  nozzles  of  very  restricted 
geometry.  One  method  is  to  reduce  the  problem  to  that  of  solving  a set  of 
ordinary  differential  equations  by  expanding  the  flow  variables  in  inverse 
power  series  in  the  radius  of  curvature  of  the  throat  (Ref.  8 and  Appendix  A). 
However,  this  method  is  limited  to  nozzles  with  a throat  wall  radius  of  curva- 
ture, R , greater  than  twice  the  throat  radius,  r (e.g.  Fig.  1).  This  can 
be  overcome  by  using  an  inverse  approach5,  in  which  the  boundary  geometry  is  not 
specified  but  is  obtained  from  the  solution  of  a prescribed  velocity  distribu- 
tion along  a suitable  reference  line.  The  velocity  is  then  modified  until  a 
streamline  of  the  flow  approximates  to  the  shape  of  the  nozzle.  This  method 
depends  on  the  ability  to  specify,  a priori,  the  velocity  distribution  which 
yields  the  desired  nozzle  boundary.  It  is  useful  for  nozzles  whose  geometry  can 
be  defined  by  the  radius  of  curvature  of  the  throat,  but  is  unlikely  to  give 
accurate  results  for  nozzles  of  complicated  geometry. 

The  method  adopted  was  that  described  by  Stow10,  which  can  be  applied  to  a 
nozzle  of  general  geometry.  It  is  based  on  the  "streamline  curvature"  method 
used  for  subsonic  flows. 

3. 1 Governing  equations 

The  equations  governing  the  steady  axisymmetric  motion  of  an  inviscid 
compressible  fluid  are,  in  the  coordinates  of  Fig.  1: 

continuity 


h (r  p u)  + h (r  p v) 


conservation  of  momentum 
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(3.3) 


Equation  (3.1)  implies  the  existence  of  a stream  function  ip  satisfying 
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By  introducing  the  speed  of  sound,  defined  by  c ■ (3P/3p)  , using  the 

relationship  v * u tan  0 , and  the  identity 
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it  is  possible  to  rewrite  equations  (3.1)  to  (3.3)  as 
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Eliminating  (3P/3x)^  between  equations  (3.5)  and  (3.6)  gives  the  radial 
equilibrium  equation: 
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+ u tan  0 


3 tan  0 
3r 


+ U-2  -t-an2  9 (3.7) 


relating  the  static  pressure  field  to  the  streamline  geometry.  The  term 
is  referred  to  as  the  "curvature’1  of  a streamline. 

♦ 


♦Subscripts  i|/  , r , x refer  to  differentiation  along  a line  of  constant  ^ , r 
or  x . 


3 tan  0 
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3 . 2 Method  of  solution 

Equation  (3.7)  is  solved  numerically  in  an  iterative  manner.  The  velocity 
field  and  the  shape  of  the  streamlines  are  approximated  from  either  a series  or 
a one -dimens ion al  solution  (see  Appendix  A).  Since  all  other  quantities  are 
known,  equation  (3.7)  then  becomes  an  ordinary  differential  equation  in  the 
static  pressure,  and  if  an  approximation  is  made  to  the  mid-radial  pressure,  it 
may  be  integrated  to  the  wall  and  axis.  The  mass  flux  is  determined  from 
equation  (3.4)  and  the  mid-radial  pressure  is  changed  until  the  calculated  mass 
flux  agrees  with  the  true  mass  flux  (see  Appendix  B) . All  the  axial  stations 
are  solved  in  this  way,  and  the  other  flow  variables  are  evaluated  from  isen- 
tropic  relations.  A new  approximation  to  the  streamline  geometry  is  made  and 
the  calculation  is  repeated  until  the  solution  converges. 


Equation  (3.7)  has  an  indeterminacy  at  a sonic  point.  If  the  values  of  the 
calculated  streamline  slope  and  curvature  were  exact  the  equation  would  have  a 
zero  on  both  sides  and  the  indeterminacy  could  be  removed.  However,  since  the 
procedure  is  iterative,  only  approximate  values  of  the  slope  and  curvature  are 
known,  and  the  equation  has  a singularity  at  such  a point.  This  problem  is 
overcome  by  using  equation  (3.5)  in  the  first  few  iterations  to  find  the  approx- 
imate locations  of  the  sonic  points,  the  term  (3P/Dx),  being  determined  from 
the  solution  in  the  previous  iteration.  Equation  (3.5)  is  then  used  close  to 
the  sonic  points  while  equation  (3.7)  is  used  elsewhere. 


In  practice  relaxed  values  of  the  calculated  slope  and  curvature  must  be 
used  in  certain  cases  to  ensure  convergence.  For  example. 
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The  relaxation  factor  depends  upon  the  nozzle  geometry  and  the  closeness  of  the 
initial  estimate  to  the  final  solution,  as  does  the  number  of  steps  required, 
convergence  being  assumed  when  the  change  in  the  flow  parameters  between 
successive  iterations  is  smaller  than  some  desired  degree  of  accuracy. 
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♦Subscripts  x and  r denote  differentiation  w.r.t  x or  r 
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particle  heat  transfer  equation 
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The  drag  and  heat  transfer  factors,  f and  g , vary  with  local  flow 

P P 

field  properties  in  accordance  with  empirical  relationships.  The  empirical 
formulae  used  are  discussed  in  Appendix  C. 

4 . 2 Basic  assumptions 

The  following  basic  assumptions  are  made  in  deriving  equations  (4.1)  to 
(4.9)  (Ref.  9): 

1)  there  are  no  mass  or  energy  losses  from  the  system; 

2)  the  gas  is  inviscid  except  for  its  interactions  with  the  condensed 
particles; 

3)  the  volume  occupied  by  the  condensed  particles  is  negligible; 

4)  the  thermal  (Brownian)  motion  of  the  particles  does  not  contribute  to 

the  pressure  of  the  system; 

5)  the  condensed  particles  do  not  interact  directly; 

6)  the  drag  and  heat  transfer  characteristics  of  an  actual  shape  and 

size  distribution  of  particles  can  be  represented  by  spherical 
particles  of  a single  size; 

7)  the  internal  temperature  of  the  particles  is  uniform; 

8)  energy  exchange  between  the  gas  and  the  particles  is  controlled  only 
by  convection; 
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9) 


the  only  forces  on  the  particles  are  viscous  drag  forces; 

10)  there  is  no  mass  transfer  between  the  gas  and  the  condensed  phase 
during  the  nozzle  expansion; 

11)  the  particles  do  not  undergo  a phase  change  in  the  region  of  the 
calculation; 

12)  the  gas  phase  is  thermally  and  calorically  perfect  and  the  flow  is 
chemically  frozen. 

The  possibility  of  eliminating  some  of  these  limitations  is  discussed  in 
Section  6. 


4 . 3 Method  of  solution 

An  initial  flow  field  is  computed  from  the  TEPGAS  program  using  a modified 
molecular  weight  (see  equation  (2.1))  and  a modified  isentropic  exponent  (see 
equation  (2.2)).  The  gas  and  particle  velocities  and  temperatures  are  then 
equated  to  the  TEPGAS  values,  and  the  gas  and  particle  densities  are  determined 
from  knowledge  of  the  mixture  density  and  particle  loading.  It  is  assumed  that 
the  static  pressure  field  and  the  locations  of  the  gas  streamlines  remain 
constant  at  the  TEPGAS  values,  and  that  the  gas  velocity  and  density  gradients 
do  not  vary  from  their  initial  values. 

The  partial  differential  equations  (4.3)  to  (4.9),  which  are  quasi-linear 
and  of  first  order,  are  rewritten  in  finite  difference  form  as  linear  algebraic 
equations,  and  used  to  derive  a set  of  algebraic  replacement  equations  in  which 
a dependent  variable  is  isolated  on  the  left  hand  side  of  the  equation  and  its 
value  determined  numerically  from  the  right. 

The  equations  are  then  solved  iteratively  by  performing  the  following 
sequence  of  operations. 

1)  By  rewriting  equations  (4.7)  and  (4.8),  the  particle  velocities  are 
calculated  from 
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The  quantities  on  the  right  of  the  equations  are  evaluated  from  the 
initial  flow  field  values. 

2)  The  particle  trajectories  are  integrated  through  the  particle  velocity 
field  via  the  equation 


tan  6 = v /u 

P P P 


and  the  new  particle  streamlines  determined.  The  outermost  particle 
streamline  constitutes  a boundary  between  the  two-phase  and  particle- 
free  regions.  Interpolation  is  carried  out  to  find  the  particle 
properties  at  the  new  particle  grid  points. 

3)  The  particle  density  is  calculated  from  equation  (4.2)  to  ensure 
continuity . 

4)  The  gas  properties  at  the  particle  grid  points  are  found  by  inter- 
polation and  the  particle  temperature  is  calculated  from  equation 
(4.9). 


5)  Steps  1 to  4 are  repeated. 

6)  The  gas  velocity  components  and  the  gas  temperature  and  density  at 
the  gas  phase  coordinate  system  grid  points  are  calculated. 


7)  The  gas  phase  total  pressure  p and  total  temperature  T are 

8o  8o 

calculated  for  each  gas  streamline  at  its  intersection  with  the 
limiting  particle  streamline.  These  quantities  remain  constant  along 
the  gas  streamlines  in  the  particle-free  region,  in  which  the  gas  is 
assumed  to  be  strictly  adiabatic,  and  this  permits  computation  of  the 
gas  velocities,  temperature  and  density  from  isentropic  relationships. 


8)  This  terminates  one  iteration  of  the  relaxation  process,  and  the 

solution  is  tested  for  convergence.  Generally,  absolute  convergence 
is  not  possible,  and  a decrease  in  the  change  in  all  variables 
between  successive  iterations  is  considered  equivalent  to  convergence. 
If  this  has  not  been  reached,  the  gas  properties  are  found  at  the 
particle  coordinate  system  grid  points  and  the  sequence  of  operations 
is  repeated. 
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RESULTS  OF  COMPUTATIONS 


The  equilibrium  particle/gas  mixture  program  (TEPGAS)  and  the  non- 
equilibrium gas/particle  program  (TGPF)  were  run  with  the  input  data  shown  in 
Tables  1 and  2.  The  nozzle  chosen  (Fig.  1)  has  a relatively  large  normalised 
throat  wall  radius  of  curvature  (>5),  but  even  so  the  results  from  TEPGAS 
(Fig.  2 and  Table  3)  show  significant  variation  in  the  flow  variables  across  the 
nozzle,  which  becomes  more  pronounced  as  the  calculation  proceeds  downstream. 

This  two-dimensional  effect  produces  the  bent  sonic  line  shown  in  Fig.  3, 

Curve  I,  which  corresponds  to  a fairly  large  transonic  region.  When  account  is 
taken  of  the  two-phase  effects,  the  particle  properties  change  very  little  from 
the  equilibrium  values  (Table  4).  There  is  a slight  increase  in  the  particle 
temperatures  over  the  TEPGAS  values,  and  the  particle  axial  velocities  are 
marginally  lower.  There  is  also  very  little  change  in  the  gas  properties  away 
from  the  nozzle  wall,  the  gas  temperatures  falling  slightly  from  the  equilibrium 
values,  and  the  axial  velocities  increasing  slightly.  However,  there  is  a very 
marked  change  in  the  gas  properties  close  to  the  nozzle  wall.  Fig.  3 compares 
the  gas  sonic  line  predicted  by  the  TEPGAS  program  (Curve  II),  and  that  given 
after  three  iterations  of  TGPF  (Curve  III)  . The  equilibrium  assumption  gives 
the  sonic  line  location  on  the  axis  quite  well,  but  provides  a very  poor  estim- 
ate near  the  wall.  This  is  because  the  particles  depart  from  the  gas  stream- 
lines, leaving  a particle-free  region  adjacent  to  the  wall.  This  departure  is 
only  slight  close  to  the  axis,  but  becomes  more  significant  where  the  curvature 
of  the  gas  streamlines  is  greater.  Figs  4 and  5 show  the  effect  of  the  isen- 
tropic  nature  of  the  gas  expansion  in  the  particle-free  region,  the  gas  tempera- 
ture and  axial  velocity  differing  by  up  to  15  per  cent  from  the  values  calculated 
with  the  equilibrium  assumption. 

The  programs  were  also  run  for  a range  of  nozzle  shapes  and  particle  sizes, 
and  for  different  particle  loadings.  A decrease  in  the  normalised  throat  wall 
radius  of  curvature  of  the  nozzle  leads  to  increased  variation  in  the  properties 
across  the  nozzle,  and  results  in  a larger  particle-free  region  close  to  the 
nozzle  wall.  The  particle-free  region  is  also  extended  if  bigger  particles  are 
introduced  (Fig.  6).  With  1 pm  diameter  particles,  the  flow  is  very  close  to 
gas/particle  equilibrium,  but  as  the  particle  size  is  increased,  the  limiting 
particle  streamline  departs  from  the  nozzle  wall  further  upstream,  changing  the 
local  stagnation  conditions  in  the  isentropic  region,  and  resulting  in  the  gas 
sonic  velocity  close  to  the  wall  being  reached  much  sooner.  This  leads  to 
greater  differences  between  the  gas  properties  predicted  with  the  equilibrium 
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assumption  and  those  calculated  taking  into  account  the  particle  velocity  and 
thermal  lags  (Fig.  7).  Changing  the  particle-to-gas  mass  ratio  alters  the 
position  of  the  gas  sonic  line.  If  the  number  of  particles  is  reduced  the  gas 
sonic  line  moves  closer  to  the  mixture  sonic  line,  but  if  it  is  increased  the 
gas  sonic  line,  in  the  two-phase  region,  is  displaced  further  downstream,  thus 
extending  the  transonic  region. 

6 CONCLUSIONS 

Two  computer  programs  to  solve  the  two-dimensional,  two-phase  flow  equations 
in  the  transonic  throat  region  of  a rocket  exhaust  nozzle  have  been  developed. 

The  first  treats  the  gas/particle  mixture  as  a heavy  perfect  gas  with  modified 
isentropic  exponent  and  molecular  weight,  and  solves  the  equations  of  isentropic 
transonic  flow  through  an  axisymmotric  nozzle.  An  initial  flow  field  estimate 
is  obtained  for  input  to  the  second  program,  which  incorporates  the  non- 
equilibrium effects  of  the  two-phase  fluid.  The  equations  governing  the  transonic, 
two-dimensional  flow  of  a gas/particle  mixture  are  expressed  as  finite  difference 
replacement  equations  and  are  solved  by  a numerical  relaxation  technique.  The 
method  is  limited  by  the  use  of  finite  difference  approximations,  and  by  the 
assumptions  listed  in  Section  4.2.  Future  work  will  be  concerned  with  lifting 
some  of  these  restrictions.  The  assumptions  that  the  particles  do  not  undergo  a 
phase  change  and  are  of  a single  size  will  be  studied  first  as  these  are  of  some 
importance  and  may  be  dealt  with  relatively  easily.  The  next,  and  most  difficult, 
step  will  be  to  remove  the  restriction  of  a frozen  gas  composition.  Equilibrium 
chemistry  (infinite  reaction  rates)  will  be  assumed  in  the  transonic  throat 
region  of  the  nozzle.  When  account  has  been  taken  of  the  chemical  reactions, 
the  possibility  of  allowing  for  particle  interactions,  and  for  the  transfer  of 
mass  between  the  gas  and  the  condensed  phase,  will  be  considered. 

The  programs  have  been  used  in  exemplary  calculations  of  the  flow  through 
an  axisymmetric  nozzle  for  a variety  of  wall  shapes,  particle  sizes  and  particle- 
to-gas  mass  ratios.  The  results  show  that  if  accurate  initial  line  data  for  the 
start  of  a supersonic  calculation  are  to  be  obtained,  the  radial  variations,  and 
the  non-equilibrium  effects  of  all  but  the  smallest  particles  (<1  urn),  must  be 
included.  When  the  modifications  discussed  above  have  been  implemented,  the 
final  stage  of  the  work  will  be  undertaken.  This  is  to  write  a computer  program 
to  solve  the  two-phase  flow  equations  in  the  supersonic  region  of  the  nozzle, 
taking  into  account  the  finite,  non-zero  rates  of  chemical  reactions. 


15 


7 


REFERENCES 


No. 

1 

2 

3 

4 

5 

6 


7 

8 

9 


Authors 

Jensen,  D.E. 
Wilson,  A.S. 

Jensen,  D.E. 
Spalding,  D.B. 
Tatchell,  D.G. 
Wilson,  A.S. 

Wilson,  A.S. 


Hoglund,  R.F. 


Hopkins,  D.F. 
Hill,  D.E. 


Kliegel,  J.R. 
Nickerson,  G.R. 


Kliegel,  J.R. 
Nickerson,  G.R. 


Kliegel,  J.R. 
Quan,  V. 


Regan,  J.F. 
Thompson,  H.D. 
Hoglund,  R.F. 


Title,  etc 

Prediction  of  rocket  exhaust  flame  properties. 
Combustion  and  Flame  1975,  £5,  43-55 

Analysis  of  recirculating  rocket  exhaust  flows. 
To  be  published 


A user's  guide  to  computer  programs  for  the 
calculation  of  conditions  in  flames  and  rocket 
nozzles . 

RPE  Tech.  Report  No.  72/10  (1972) 

Recent  advances  in  gas-particle  nozzle  flows. 

ARS  Journal,  May  1962,  662-671 

Effect  of  small  radius  of  curvature  on  transonic 
flow  in  ax i symmetric  nozzles. 

AIAA  Journal,  1966,  4,  (8),  1337-1343 

Flow  of  gas-particle  mixtures  in  axially  symmetric 
nozzles . 

Progress  in  Astronautics  and  Rocketry,  v.6. 
Detonation  and  two-phase  flow;  edited  by 
S.S.  Penner  and  F.A.  Williams. 

Academic  Press,  New  York,  1962,  pp.  173-194 

Axisyrametric  two-phase  perfect  gas  performance 
program. 

TRW  Systems  Group,  Redondo  Beach,  California, 

Rept.  02874-6006-R000,  Vols  I & II  (1967) 

Convergent-divergent  nozzle  flows. 

TRW  Systems  Group,  Redondo  Beach,  California, 

Rept.  02874-6002-R000  (1966) 

Two-dimensional  analysis  of  transonic  gas-particle 
flows  in  axisymmetric  nozzles. 

J.  Spacecraft  and  Rockets,  1971,  8,  (4),  346-351 


16 


No. 

Authors 

Title,  etc 

10 

Stow,  P. 

The  solution  of  isentropic  transonic  flows. 
J.  Inst.  Maths  Applies,  1972,  9,  35-46 

11 

Carlson,  D.J. 
Hoglund,  R.F. 

Particle  drag  and  heat  transfer  in  rocket  nozzles. 
AIAA  Journal,  1964,  2,  (11),  1980-1984 

12 

Crowe,  C.T. 

Drag  coefficient  of  particles  in  a rocket  nozzle. 
AIAA  Journal,  1967,  5,  (5),  1021-1022 

13 

Crowe,  C.T. 

Inaccuracy  of  nozzle  performance  predictions 
resulting  from  the  use  of  an  invalid  drag  law. 

J.  Spacecraft  and  Rockets,  1970,  (12),  1491-1492 

14 

Kliegel,  J.R. 

Gas  particle  nozzle  flows.  Ninth  International 

Symposium  on  Combustion. 

Academic  Press,  New  York,  1963,  pp.  811-827 


17 


APPENDIX  A 


Approximate  solution  of  isentropic  transonic  flow  by  series  solution  method 

This  method,  developed  by  Kliegel  and  Quan^®,  depends  upon  the  nozzle  wall 
shape  being  expressed  in  terms  of  R , the  nozzle  throat  wall  radius  of  curva- 
ture, normalised  with  respect  of  the  nozzle  throat  radius  rfc  . Thus,  an  R 
is  chosen  such  that  the  nozzle  wall  may  be  approximated  by  the  hyperbola 


1 + z' 


where  z 


Ji  r 


— and  y ■ — . 

t rt 


The  axial  and  radial  velocity  components,  u and  v , may  be  expressed 
in  inverse  powers  of  R by  the  expansions 


ux  (y , z)  u2  (y , z) 

u - uo  (y,z)  + + 5 + 

R 


V,  (y,z)  V (y,z) 

v (y,z)  + + + 

o '•7*  R r2 


It  may  be  shown  that 


where 


uq  (y,z)  - aQ  (z) 


vQ  (y,z)  - aL  (z).y 


a dy 
o v 

*1  y dz 
'w 


(1  - ./) 


aa  . _ \ a dv 

+ 2 1 -1^-a  2 -o  « - 

dz  Y+loJydz 


with  boundary  conditions 
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since  at  the  throat 


aQ  (0)  «*  1 


ax  (0)  = 0 , 


dy, 

dz 


w 


Similarly,  the  first  order  velocity  components  are  defined  by 


ux  (y,z)  = bQ  (z)  + b2  (z)  y" 


Vj  (y,z)  = bx  (z)  y + b3  (z)  y 


where  bQ  - b^  are  found  by  solving  four  ordinary  differential  equations 
similar  to  those  relating  aQ  and  a^  . These  equations  are  algebraic  at  the 
throat  and  can  be  solved  directly  to  find  the  throat  boundary  conditions.  The 
second  order  velocity  components  are  defined  by 


u2  (y.z)  = (z)  + C2  (z)  y2  + C4  (z)  y4 


v2  (y,z)  = Cj  (z)  y + C3  (z)  y3  + C5  (z)  y5 


where  Cq  - are  related  by  six  ordinary  differential  equations  which  can 
also  be  solved  algebraically  to  determine  the  throat  boundary  conditions. 

For  a nozzle  whose  shape  may  be  approximated  by  a hyperbola  with  a normal- 
ised radius  of  curvature  greater  than  three,  the  second  order  series  solution 
yields  a very  good  first  approximation.  However,  for  nozzles  of  sharp  wall 
curvature  the  estimates  are  not  so  reasonable,  particularly  if  R is  less  than 
two.  For  a normalised  throat  wall  radius  of  curvature  of  less  than  one,  the 
second  order  solution  predicts  that  the  throat  axis  velocity  is  supersonic, 
which  is  physically  impossible.  The  third  and  higher  order  equations  may  be 
obtained,  but  the  mathematics  involved  becomes  extremely  cumbersome.  Thus,  for 
nozzles  where  the  second  order  solution  is  not  feasible,  a one-dimensional 
approximation  is  used. 
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APPENDIX  B 


Determination  of  the  maaa  flu?. 

The  mass  flux  m is  calculated  at  the  throat  in  each  iteration  by  choosing 
the  mid-radial  static  pressure  which  maximises  the  mass  flux  through  the  throat. 
This  is  equivalent,  in  one-dimensional  flow,  to  finding  the  point  where  the  Mach 
number  becomes  unity. 

At  the  throat 


m 


(r  p u) 


x const 


dr 


(B.l) 


By  writing 


P ■ PQ  (P/Pq)^Y  ■ pq  where  z 


P/P 

o 


and 


u ■ Me  cos  0 


•I 


2y  R T 
c 

y - l 


cos  e 


1/2 


equation  (B.l)  may  be  converted  into 


Cl 


m ■ 2ir  K { r cos  0 zJ 
o 


i/r  f t , Jd 


n 


dr 


where 


2T  P0  1 

Y - 1 


For  m to  be  a maximum,  z is  such  that 


37  ■ 2*  * J‘  [«2/t  - «1r]'1/2[2'T  *J/T‘1  - *T  ll/l]  ( or  - 0 
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When  the  required  mass  flux  has  been  found,  the  other  axial  stations  are 
solved  by  determining  the  mid-radial  pressure  which  satisfies 


f (z) 


(B.2) 


Since  equation  (B.2)  has  a turning  point  at  the  throat,  at  axial  stations 
upstream  and  downstream  of  the  throat  two  possible  solutions  exist.  These 
correspond  to  the  subsonic  and  supersonic  solutions  of  one-dimensional  flow, 
although  two-dimensional  flow  can  be  mixed  in  the  radial  direction.  The  problem 
of  which  solution  to  choose  arises  only  at  stations  downstream  of  the  throat, 
since  upstream  the  equivalent  of  the  one-dimensional  subsonic  solution  occurs. 
Physically  there  is  continuous  change  in  the  mid-stream  Mach  number  for  isen- 
tropic  flows  which  serves  as  a basis  for  choosing  the  solution  downstream  of 
the  throat. 
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APPENDIX  C 


Particle  drag  and  heat  transfer  factors 


Definitions 

The  empirical  relationships  used  to  calculate  the  particle  drag  and  heat 
transfer  factors  are  the  same  as  those  used  by  the  TRW  program ^ and  by  Regan^. 
They  are 

fp  * *D  (1  * 2.52  C)/(l  + 3.78  C) 
and 

gp  - Kp/(1  + 3.42  Kp  c/y1/2  Pr) 


where 


V / P r (R  T ) 
g 8 P g 


1/2 


and 


*0  “ Stokes 


where  Cp  is  the  particle  drag  coefficient. 

The  gas  viscosity  is  a function  of  temperature 


- u (T„/T  ) 

8 8o  8 8o 


0.67 


and  the  Prandtl  number  is  computed  from  Euken's  relationship 


Pr  ■ 4y/(9y  - 5) 


Drag  coefficient 

Because  the  condensed  particles  are  assumed  to  be  spherical  the  drag 
coefficient  for  a sphere  is  used  to  estimate  the  particle  drag  coefficient. 
Experimental  work  to  establish  the  magnitude  of  the  drag  forces  on  spheres  has 
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been  going  on  since  Newton  began  his  experiments  in  1710,  but  until  recently  there 
were  no  experimental  data  to  give  the  drag  coefficient  of  a sphere  in  the  flow 
regime  of  interest.  In  1850  Stokes  solved  the  Navier-Stokes  equation  by 
neglecting  the  inertial  terms  (Re  -*•  0)  and  obtained  the  simple  drag  law 

<°D>Stok..  * 24  /R'  <C'l> 

where  Re  is  the  particle  Reynolds  number  based  on  the  speed  of  the  particle 
relative  to  the  gas  and  is  defined  by 

Re  = 2 r (W  - W)  p /u  . (C.2) 

P g P g g 

* However,  it  has  been  shown4  that  the  Stokes  drag  law  has  only  limited 
applicability  to  gas/particle  nozzle  flows  since  it  is  restricted  to  continuum, 
incompressible  flow  and  particle  Reynolds  numbers  less  than  1.  The  flow  regimes 
encountered  by  the  micron-sized  particles  extend  from  continuum  to  free-molecule 
flow,  and  particle  Reynolds  numbers  of  up  to  100  and  particle  Mach  numbers 
exceeding  1.  It  was  therefore  necessary  to  derive  expressions  which  agreed  with 
the  available  data,  conformed  with  theoretically  predicted  trends,  and  provided 
reasonably  smooth  variations  with  Mach  and  Reynolds  numbers. 

More  recently,  drag  coefficient  data  have  become  available,  and  Crowe13 
compares  the  equations  devised  by  Kliegel14,  Carlson11  and  Crowe12  with  an 
empirical  expression  which  he  developed  to  fit  the  new  data.  He  concludes  that 
the  differences  between  nozzle  performance  predictions  made  with  each  devised 
drag  law,  and  those  made  using  the  law  validated  by  experiment,  are  negligible. 
Thus,  since  the  choice  of  drag  law  has  little  effect  on  the  predicted  results, 
it  was  decided  to  use  the  relationship  derived  by  Carlson11  as  this  is  the 
simplest  expression.  It  is  based  on  the  Stokes  drag  law,  equation  (C.l),  and 
includes  terms  to  correct  for  rarefaction,  inertial  effects  and  compressibility 
effects: 


(1  ♦ 0.15  Re0,687)  [ 1 ♦ exp  (-0 .427/M4 *63  - 3/Re0,88)  ~ 

1 ♦ (M/Re)  Q 3.82  + 1.28  exp  (-1.25  Re/M)  ] 

where  Re  is  the  particle  Reynolds  number  and  M is  the  particle  Mach  number 
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Table  1 


Input  data  for  TEPGAS* 

Throat  radius,  r 

0.4328  cm 

Radius  of  curvature  of  circular  arc  section,  R 

2.54  cm 

Distance  from  nozzle  throat  to  start  of  conic  section 

0.6574  cm 

Angle  of  conic  section 

15° 

Modified  isentropic  exponent,  7 

1.157 

, -2 

Chamber  pressure,  Po 

3.507  E 06  Nm 

Chamber  temperature,  Tq 

2472  K 

Modified  gas  constant,  R/m 

239.0  Nm  K_1  kg"1 

No.  of  streamlines 

16 

No.  of  axial  stations 

40 

Axial  step  length 

0.04  cm 

Table  2 


Input  data  for  TGPF* 

Gas  isentropic  exponent. 

Yg 

1.232 

Gas  reference  viscosity, 

P 

8.0  E-05  kg  m 1 s 1 

Gas  reference  temperature 

go 

♦ To 

2500  K 

Gas  specific  heat,  Cp^ 

80 

1500  J kg"1  K_1 

Particle  radius,  r 

1 .0  E-06  m 

P 

Particle  density,  m 

D 

3.97  E 03  kg  m"3 

Particle  specific  heat. 

Cpp 

1355  J kg"1  K-1 

Particle-to-gas  mass  flow 

ratio,  <(> 

3/7 
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Nomenclature 


c 


C 

P 

F, 

G, 


D 


f 

P 


m 

m 

m 

P 

m 

M 

P 

Pr 

r 

r 

P 

R 


Re 

T 

u 

v 

W 

x 


Y 

Y 

e 


u 

p 

p 


g 

P 


4> 

<t> 


speed  of  sound 
particle  drag  coefficient 
specific  heat  at  constant  pressure 
particle  drag  factors 
particle  heat  transfer  factors 

Stokes 

molecular  weight 

molecular  weight  of  gas  particle  mixture 

density  per  unit  volume  of  particle 

mass  flow  rate 

gas  Mach  number 

pressure 

Prandtl  number 

radial  coordinate  measured  from  nozzle  axis 
particle  radius 

gas  constant;  nozzle  throat  wall  radius  of  curvature,  normalised  with 
respect  to  throat  radius  r 

Reynolds  number 

temperature 

axial  velocity  component 
radial  velocity  component 
speed 

axial  coordinate  measured  from  nozzle  throat 
isentropic  exponent 

isentropic  exponent  of  gas/particle  mixture 
streamline  angle  with  respect  to  nozzle  axis 
viscosity  coefficient 
gas  density 

particle  density  in  the  gas  (based  on  the  gas  volume) 
stream  function 

particle-to-gas  mass  flow  rate 


Subscripts 


g gas  property 

o total  condition 

p particle  property 

w nozzle  wall  condition 

t throat  condition 
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MACH  NUMBER  NOZZLE  RADIUS,  NORMALISED  BY 


MEMO  685 
FIG.  3 


GAS  AXIAL  VELOCITY, m sec'1  GAS  T EMPERATURE,K 


MEMO  685 
FIG.  U & 5 


FIG. 4 COMPARISON  OF  GAS  TEMPERATURES  ALONG  NOZZLE 
WALL  FOR  2(jm  PARTICLES. 


FIG. 5 COMPARISON  OF  GAS  AXIAL  VELOCITIES  ALONG 
NOZZLE  WALL  FOR  2pm  PARTICLES. 


DISTANCE  from  NOZZLE  WALL  of 

GAS  AXIAL  VELOCITY  m sec"1  LIMITING  PARTICLE  STREAMLINE , cm  kIC  03 
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FIG.  6 «7 


FIG.6  DISPLACEMENT  FROM  THE  NOZZLE  WALL  OF  THE 
LIMITING  PARTICLE  STREAM  LINE  , FOR  DIFFERENT 
SIZE  PARTICLES. 


FIG. 7 COMPARISON  OF  GAS  AXIAL  VELOCITIES  ALONG  THE 
NOZZLE  WALL  FOR  DIFFERENT  SIZE  PARTICLES. 
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Abet  ran  This  memorandum  describes  progress  in  the  development  of  a set  of  computer 
programs  which  solve  the  equations  of  gas/particle  flow  through  an  axisyiumetric 
nozzle.  Two  programs  have  been  written  to  solve  the  equations  in  the  transonic  throat 
region  of  the  nozzle.  The  first  treats  the  two-phase  fluid  as  a heavy  perfect  gas 
with  modified  isentropic  exponent  and  molecular  weight,  and  solves  the  equations  of  | 
isentropic  transonic  flow.  An  initial  flow  field  configuration  is  obtained  for  input  ' 
to  the  second  program,  which  incorporates  non-equilibrium  effects  of  the  gas/particle  ! 
mixture . 
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The  programs  provide  accurate  initial  line  data  for  starting  a supersonic  calculation, 
and  will  be  used  as  input  to  a third  program,  soon  to  be  written,  for  solving  the  two-| 
phase  flow  equations  in  the  supersonic  region  of  the  nozzle.  * | 
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